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ABSTRACT 



A transient overpower (TOP) accident in a Liquid Metal 
Fast Breeder Reactor (LMFBR) is considered. The analysis 
is formulated to model the dynamic response of the reactor 
fuel subassembly during the initial period of the postulated 
overpower transient. An equivalent cylindrical cell is used 
to model the fuel subassembly. The governing neutronic and 
heat transport equations for each region (fuel, clad, and 
coolant) of the equivalent cylindrical cell are developed. 
Nuclear Doppler broadening feedback is included in the dynam- 
ic model making the coupled equations non-linear. The re- 
sulting non-linear partial different i a 1 field equations are 
transformed into a system of ordinary differential equations 
by the finite element method. An isoparametric, quadratic, 
rectangular element is used for the discretization of the 
spatial domain. When using the finite element method, large 
system matrices may result. To facilitate solution of these 
large systems, an optimum compacting scheme is utilized. 

The implicit Gear's method is used for the solution of the 
system of ordinary differential equations. The results for 
a sample problem are presented. 
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I. INTRODUCTION 



As the world's fossil fuel resources are depleted, more 
emphasis is being placed on the breeder reactor as a poten- 
tial means of solving the coming energy crisis. While the 
development of new energy sources i s being pushed, equal ef- 
fort is being given to the maintenance of an environmental- 
ly clean world. To this end, the safety of breeder reactors 
is receiving a considerable amount of attention before assum- 
ing that the breeder reactor is the answer to the energy 
problem. 

The Liquid Metal Fast Breeder Reactor (LMFBR) appears to 
be one of the most promising breeder reactors. Most engineers 
will concede there is little probability of a nuclear explo- 
sion occurring in the operation of a nuclear reactor. Of 
major concern to engineers is the loss of coolant accident 
(LOCA) and the transient overpower accident (TOP). The pres- 
ent analysis is concerned with a TOP accident in a LMFBR. 

The analysis is formulated to model the dynamic response of 
the reactor fuel subassembly during the initial period of the 
postulated overpower transient. The primary consideration is 
given to the early response of this fuel subassembly to 
various conditions of disturbances. The phenomenon which oc- 
curs after core disassembly (i.e., clad melting) is not the 
concern of this analysis. Only the time prior to clad melting 
is being considered. 
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No consideration is given here as to how the overpower 
transient occurs or to why the safety features of the reac- 
tor did not operate properly. It is postulated that the 
accident has occurred. In this analysis, the TOP accident 
is created by either a step increase in reactivity, a ramp 
increase in reactivity, or a combination of both. 

An inherent safety feature of most reactors, nuclear 
Doppler broadening feedback, is included in the dynamic model 
of the fuel subassembly. The Doppler feedback acts to reduce 
the effect of the excursion. Consideration of this feedback 
creates a non-linear system model which is described by a 
non-linear, i ni ti al - boundary- val ue problem. 

The conventional method of solution uses the standard 
point kinetics formulation. Recent studies have pointed out 
a non-negl i gi bl e error in this model [1], particularly with 
asymmetric disturbances [2], or space-dependent feedback [3], 
In Ref. [4], a somewhat novel approach of using the finite 
element method (FEM) for the space-time dependent solution 
of the reactor dynamics problem was demonstrated. The FEM 
is effective in handling these asymmetric disturbances and 
space-dependent feedbacks. Therefore, the finite element 
method was used so that the spatial effects on the postu- 
lated problem may be studied further. 

The purpose of this work was to demonstrate further the 
applicability of the FEM to the non-linear reactor dynamics 
problem as well as to i nvesti gate the dynamic response of the 
reactor fuel subassembly. The analysis required a novel 
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approach to handle the gap conductances present at the in- 
terfaces of the equivalent cell model of the fuel subassem- 
bly; this will be clarified in the analysis. 
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II . DESCRIPTION OF PROBLEM 



A. PHYSICAL SYSTEM 

The typical Liquid Metal Fast Breeder Reactor (LMFBR) 
core consists of many hexagonal modules, each containing 
several hundred fuel pins. For this analysis, an equivalent 
cylindrical cell is used to model the fuel subassembly; see 
Figure 1. The use of equivalent cells as models for larger 
systems has been common practice in nuclear analysis (i.e., 
the well known Wigner-Sietz method). In using an equivalent 
cell, the actual shape of the reactor core is not important, 
and the analysi s is applicable to any reactor which has the 
same equivalent cell. 

The equivalent cell considered in this analysis, Figure 1, 
is fueled with enriched uranium dioxide, has a stainless 
steel cladding, and has liquid sodium for a coolant. The 
dimensions used are 

a = 0.254 cm, 
b = 0.292 cm , 
c = 0.365 cm, 

and H = 33.0 cm. 

The gap between the fuel and cladding is very small and, in 
fact, may be nonexistent as in bonded fuels. The dimension 
of this gap has been assumed negligible. The height, H, of 
the fuel rod is shorter than many proposed systems (Fast Flux 
Testing Facility and Clinch River Breeder Reactor). However, 
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to facilitate the numerical solution a smaller rod was used. 
The dynamic behavior prior to fuel pin failure for this sys- 
tem should be similar to the behavior of larger systems. 

The treatment of this problem in three dimensions would 
be prohibitive in computer usage. Therefore, azimuthal sym- 
metry is assumed, and the problem becomes a two-dimensional 
cylindrical (r,z) problem. 

B. SYSTEM MODEL 

The analysis considers the monoenergeti c neutron diffu- 
sion approximation to model the transient neutron transport 
problem. A simple conduction-convection heat transfer model 
is used for the energy transport problem. The temperatures 
in the model are directly coupled to the neutron flux through 
the nuclear heat generation within the fuel. The neutron 
population is in turn coupled to the temperature through any 
of a number of reactivity feedback mechanisms. The nuclear 
Doppler effect is perhaps the most important of these mecha- 
nisms since it provides a negative temperature coefficient 
which increases the inherent stability of the reactor. Prior 
to core disassembly and fuel melting, the nuclear Doppler 
effect is the most dominant feedback and is, therefore, the 
only feedback mechanism considered in this analysis. 

For irradiated, mixed-oxide fuels, a phenomenon of fuel 
restructuri ng has been commonly observed. This res tructuri ng , 
essentially a change of phase of the fuel material, presents 
a unique heat transfer problem particularly during transient 
conditions. The problem has not been fully characterized and 
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is beyond the scope of this analysis. The fuel was, there- 
fore, assumed to be a homogeneous mixture of enriched uranium 
dioxide. 

At the fue 1 -cl addi ng interface, there exists a gap which 
produces a thermal resistance. This thermal resistance is 
one of the most significant deterrents to the energy transfer 
to the coolant. The interface may be in physical contact or 
an actual gap may exist. The prediction of the thermal resis- 
tance is extremely complicated and must take into considera- 
tion many parameters: initial dimensions, type of bond, fill 
gas composition, fuel restructuring, fuel swelling, prior fuel 
life cycle, to mention a few. Reference [5] documents a com- 
puter program which attempts to predict the gap conductance, 

H__ . This treats the thermal resistance at the interface in 
gap 

the same manner as a convection heat transfer coefficient when 
considering convection heat transfer. Since it is not the 
objective of this analysis to predict the gap coefficient, a 
representative set of values for gap coefficient, as given in 
Ref [6], is used in this analysis. These values are assumed 
to remain static during the transient. The gap conductance 
profile actually varies with time and will have an effect upon 
the transient, as noted in Ref. [7]. The prediction of this 
variance was not considered important for this analysis; 
therefore, the static assumption was made. 

An average convection heat transfer coefficient was used 
to determine the heat transfer from the cladding to the cool- 
ant. The value used was determined from an empirical formula 
given in Ref. [5] and repeated in Appendix C. 
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In this work, consideration is given to step and ramp in- 
creases in reactivity, although any reactivity transient may 
easily be considered. The step and ramp increases in reactiv- 
ity probably represent the most realistic physical reactivity 
inputs in a reactor. Once the reactivity has been inserted, 
the transient overpower excursion begins. Unless the Doppler 
feedback can override the inserted reacti vi ty, the excursion 
will continue until there is physical core disassembly. 

C. NUMERICAL SOLUTION 

The system of equations which models the proposed problem 
is a non-linear, i ni ti al -boundary-val ue problem. The conven- 
tional method of solution of the reactor dynamics is the point 
kinetics formulation. It was pointed out in Refs. [1], [2], 
[3], and [4], that there is a non-negligible error in this 
model, particularly under conditions of asymmetric disturbances 
or space-dependent feedback. Reference [4] demonstrates the 
somewhat novel approach of using the finite element method 
(FEM) to solve the space-time dependent reactor dynamics prob- 
lem. As shown in Ref. [4], the FEM is quite effective in 
handling localized perturbations and space-dependent feedback. 
In this work, only uniform disturbances were considered; how- 
ever, the feedback model was space-dependent. Therefore, the 
finite element method is used to solve the non-linear, coupled, 
space-time dependent neutronic and heat transport field equa- 
tions. The solution technique results in a large computer 
storage requirement; therefore, an optimum compact storage 
scheme, Ref. [8], is utilized for storage of the discretized 
matri ces . 
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Once the domai n has been discretized by the FEM, the solu- 
tion of the resulting ordinary differential equations was to 
be accomplished by the implicit Gear's method. Ref. [9]. 
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III. MODEL DEVELOPMENT 



A. NEUTRONIC ANALYSIS 

In this section the governing field equations for the 
neutron population (flux) for each of the three regions 
(fuel, clad, and coolant) of the domain will be formulated. 
The monoenergeti c, diffusion theory will be used. 

Consider an arbitrary volume of material within a reac- 
tor. Applying the condition of conservation to the mono- 



energetic neutrons leads 


to the neutron equation 


conti nui ty [10] 






3n ( r , t ) 
3 1 


= S(r,t) 


- Z a ( r ) $ ( r , t ) - div J(r, 


whe re 






r_ - 


spatial 


point 


t - 


time 




n(r,t) - 


neutron 


d e n s i ty 


S(r,t) - 


neutron 


producti on 


x a (r) - 


neutron 


absorption cross section 


<J>(r,t) - 


neutron 


flux 


J(r,t) - 


neutron 


current 



The left-hand side of equation (1) represents the time rate 
of change of the neutron density which is related to flux, 

<f> > by 

n (r,t) = 1 <*>(r,t) (2) 
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where 



v - neutron velocity 

On the right-hand side of equation (1), the first term is 
neutron production, the second term is a neutron loss through 
absorption, and the third term is a neutron loss through leak- 
age from the control volume. 

Using equation (2) and applying Fick's Law to the equation 
of continuity results in the classical neutron diffusion equa- 
tion 

V«(D(r)V4>(r,t)) -Z fl ( r ) <j> ( r , t ) +S ( r , t ) (3) 

where D(_r) - neutron diffusion coefficient 

The vector notation used here is intended to include only two 
dimensions (r,z) since azimuthal symmetry has been assumed. 

Equation (3) is applicable to each of the three regions 
of the equivalent cylindrical cell. The subscripts F (fuel), 
c (cladding), and co (coolant), will be used to denote these 
regions. 

1 . Fuel Region 

In applying equation (3) to the fuel region the 
material properties of the fuel must be used. Within the 
fuel the neutron souce term is due to the nuclear fission 
process. During fission, neutrons are released as both 
prompt neutrons and delayed neutrons so that 

S(r,t) = S p (r,t) + S D (r,t) (4) 
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where 



Sp(r.,t) - prompt neutron source 
S D (r.,t) - delayed neutron source 

The neutron sources are commonly represented as [10] 



S p ~ Z aF ^F ^ ~ ^ 



( 5 ) 



and 



'D C i X , 



( 6 ) 



where 



k oo " infinite multiplication factor 

6 - fraction of fission neutrons which appear 

as delayed neutrons 

n - number of delayed neutron groups 

r 

C.j - concentration of delayed neutron precursors 
i n the i th group 

- decay constant of the delayed neutrons 



The space and time variables, _r and t, will be dropped except 
where needed for clarification. 

The concentration of delayed neutron precursors, , 
is given by the following first order partial differential 
equation [10] 



9 C . 

Tt 1 = 6 i k «> S aF^F " X i C i 



(7) 



where 



B.j - fraction of delayed neutrons which appear 
as delayed neutrons in the i th group 
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The solution of equation (7) is 



C. = 6 i ^ VM t-t ' ) kJr,t')Z aF <l> F dt , + C i e ’ Xit (8) 
o 

C? - concentration of delayed neutron precursors 
1 of the i th group at time zero 

Reference [10] develops an expression for the initial concen- 
tration of delayed neutrons 

C i = B i k °. Z aF *° /X 1 (9 > 



where 



- initial finite multiplication factor 
4> 0 - initial neutron flux 

Combining equations (3), (4), (5), (6), (8), and (9), yields 
the governing equation for the fuel region 

V • ( DpV<j> F ) + 2 aF <J) F [k oo (l-6) -1] 

n t 

+ Z^ i C6 i f e ~^i (t-t ' ) k oo (_r,t' )4> F 2 p dt 1 
i = l o 

_ \ f i 3^ p 

'•iklWV 'i s 7r <’°> 

2 . Cladding Region 

The cladding separates the fuel and coolant and con- 
tains the fuel and the fission by-products. Equation (3) 
governs the neutron flux in the clad. Since the clad contains 
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no fissile material, the neutron source term is zero. The 
field equation is, then, 

'•("A) Vc-ir <’i) 

3 . Coolant Region 

The annular region around the cladding of the equiva- 
lent cell contains the coolant. As in the clad, there is no 
fissile material in the coolant and, therefore, no neutron 
source term. The equation governing the neutron flux in the 
coolant is 

1 3 *co 

V • ( D Vcb ) - 2 d> = - - . co (12) 

' co v co ; aco^co v 3t . ' 

Equations (10), (11), and (12) are the one- vel oci ty , 
diffusion approximation used to model the neutron transport 
probl em. 

4 . Infinite Multiplication Factor 

The infinite multiplication factor, k^, may be ex- 
pressed as the infinite multiplication factor at time zero 
(start of transient), k^, plus the postulated reactivity in- 
sertion (such as a step or a ramp), p, minus the change in 
the Doppler reactivity feedback, Ap^. Other feedback mechan- 
isms are normally not as significant as the Doppler broaden- 
ing feedback prior to fuel melting and have been neglected 
in this analysis. Therefore, 

k oo = K + P - A Pp (13) 
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For a fast reactor, k° may be approximated as 



k 



o 

00 




( 14 ) 



where 



v - average number of neutrons released 
per fission 

- fission cross section of the fuel 
I a p - absorption cross section of the fuel 

The nuclear Doppler effect is a very important safety 
feature in a nuclear reactor. Nuclei in an atom are in con- 
tinual motion due to their own thermal energy. As a result 
of this motion, even when monoenergetic neutrons interact 
with the atom, there appears to be a spread in the energy of 
the neutron - the Doppler effect. It can be shown that the 
cross section of a resonance becomes less in magnitude and 
wider as the motion of the nuclei increases [10], As the 
temperature increases, the motion increases, and the shape 
of a resonance cross section broadens. This broadening in- 
creases the average cross section, thus, providing a negative 
temperature coefficient. This effect is shown in Figure 2. 

It is this nuclear Doppler broadening effect which provides 
one of the few inherent, reliable, negative reactivity feed- 
backs which slows an overpower transient and possibly stops 
a mild overpower transient. 
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CROSS SECTION [barnt} 




Figure 2. Doppler Broadening of a Resonance Peak 
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The Doppler reactivity change with respect to fuel 



temperature changes should be written as [6] 



d P D 

dT 



aT" 3/2 + bT" 1 + cT m_1 



( 15 ) 



where a, b, and c are parameters determined from experi- 
mental work and m is an integer. However, as noted in 

Ref. [6], a substantial amount of work has shown that 
d P D 

T is very nearly constant over the temperature range 

under cons i derat ion. Therefore, the coefficients a and 
c have been set equal to zero, and b is defined as 



dp, 

K D = T dT 



( 16 ) 



The constant. Kg, is commonly called the Doppler constant. 
Solving equation (16) for pg yields 



Pg=blnTp+K 



( 17 ) 



where K is an integration which may be obtained from 
initial conditions. 



K = pj - b In T° 



( 18 ) 



where 



Pg - Doppler effect at time zero 
Tp - fuel temperature at time zero 
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Substituting for K in equation (17) will give 



P D - p° = Ap D = b 1 n (Tp/Tp) (19) 

The infinite multiplication factor now becomes 

= k° + p - b 1 n ( Tp/Tp ) (20) 



The effect of delayed neutrons is small compared to 
the prompt neutron effect; therefore, it may be assumed that 
the Doppler effect on delayed neutrons is insignificant to 
the overall problem. The of equation (8) is, then, 

assumed to be k^ . With this assumption and equation (20), 
the neutron diffusion equation in the fuel, equation (10), 
may be rewritten as 

V • ( DpV<f> F ) + Z aF 4>p[k°(l-8)-l+p(l-6)-b(l-6)ln(Tp/T°)] 



n 

+ E 

i =1 



X i (6 



i Z aF y' t e- X i (t - t,) [r + P ]<})pdt' + C?e- X i t }= 1 
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Tt 



( 21 ) 



To facilitate the present analysis, the number of 
delayed neutron groups is taken as one averaged group. So 
that, A ■ and 6 ■ become X and 6 , respectively. This approxi- 
mation should have little effect on the problem under con- 
s i derat i on . 

5 . Boundary Conditions 

The neutron diffusion problem involves solution of 
the partial differential equations (11), (12), and (21), with 
the following boundary, interface, and initial conditions: 
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1) 



0 



9<f>r 

— (o.z.t) = 

2) 4> F (a,z,t) = 

8 (j) p 

3) D F yp- (a,z,t) = 

4) <f> c (b,z,t) = 

9<f> 

5 ) D c T F~ ( b,z,t ) 

3<j> 

6) ^(c.z.t) = 

7) 4» F ( r* , +_ j,t) = 

8) ^ p ( L > 0 ) = 

9) <f> c (r,0) = 

10) <f> co (r,0) = 



4> c ( a , z , t) 

9<f> 

D c Tr“ (a- 2 - 1 ) 

* co (b,z,t) 

S^. 

D co 3F^ f 5 - 2 ’ 1 ) 

0 

<f> c (r,+ 7 »t) = 4> co (r,+ f.t) 

<t>° F (I) 

*c (r) 

*co ( ^ 



0 



Boundary condition 1) results from the assumed azimuthal sym- 
metry. Interface conditions 2), 3), 4), and 5) are continu- 
ity conditions of the flux. Boundary condition 6) results 
from the use of an equivalent cell and basically indicates 
there is an equal number of neutrons transferred in and out 
of the cell at the outer boundary. This should be valid un- 
less the cell is located near the outer edge of the reactor. 
Boundary condition 7) is an assumption that the flux is zero 
at the axial boundaries of the cell. Initial conditions 8), 
9), and 10) are the assumed initial distributions of the 
neutron flux. 
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B. HEAT TRANSFER ANALYSIS 

In this section, the principle of conservation of energy 
will be used to formulate the governing field equations for 
the heat transport in each of the three regions. A simple 
heat conduction model with convection heat transfer to the 
coolant is used to model the heat transport problem. A gap 
conductance model is used to describe the heat transport 
across the gap at the fuel-clad interface. 

1 . Fuel Region 

Conservation of energy within the fuel region yields 
the unsteady heat conduction equation with a generation term 

9Tp 

V* (k F (r)VT F (r,t) )+q(r,t)=p F (r)C pF (r) (r,t) (22) 



whe re 



kp(r) 
T p (r ,t) 
q(r,t) 

P F (r) 

c pfW 



thermal conductivity of the fuel 
fuel temperature 

nuclear energy generation per unit 
vo 1 ume 

fuel density 
fuel specific heat 



As in the neutronic analysis, the vector notation is intended 
to include only two dimensions, (r,z). The r_ and t will be 
dropped except where needed for clarification. 

The nuclear generation term may be expressed as 

q - e Zf F 4>p (23) 
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where 



e - nuclear energy released per fission 



As can be seen, it is through the nuclear generation term 
that the temperature is directly coupled to the neutron 
flux. This coupling and the temperature dependent Doppler 
reactivity feedback combine to make the coupled problem non- 
linear. 

Substituting equation (23) into equation (22) yields 
the governing thermal equation for the fuel 

V* ( kpVTp) + e Z fF <f> F = Pp C pp jjr- (24) 



2 . Cladding Region 

Conservation of energy within the clad will yield 
the heat conduction equation. With nuclear generation, a 
relatively small amount (-5%) of the energy will be released 
in the cladding and the coolant. However, in this analysis 
it is assumed that the total energy release is in the fuel 
region. This should not create any significant error. Us- 
ing this assumption, the unsteady heat conduction equation 
for the cladding becomes 



V-(k c VT c ) 



- 

P C pc 3t 



(25) 



3 . Coolant Region 

Once again, conservation of energy will lead to 
the heat conduction equation plus an additional term which 
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takes into consideration the coolant flow. The governing 
equation is 




co 



(26) 



where 



V 



co 



coolant flow velocity 



Equations (24), (25), and (26) are the governing 
equations used to model the energy transport problem. 

4 . Interface Conditions 

The interface between the fuel and the cladding may 
be an actual gap with a finite distance, or the surfaces 
may be in intermittent contact on a microscopic scale. To 
model the heat transfer across this interface, a gap heat 
transfer coefficient is introduced. The gap coefficient 
must take into consideration many items (e.g., radiation 
heat transfer across the gap, heat transfer by solid-to-solid 
contact, heat conduction across a gas filled gap). The pre- 
diction of this gap coefficient is extremely complicated and 
beyond the scope of this work. In Ref. [6], a set of values 

for H is given and the axial variation of H, „ in this 
gap s gap 

analysis is approximately the same. A cosine curve has been 
fitted to the sample data to determine the gap coefficient, 
see Figure 3. The heat flux across the fuel-clad interface 
is, then , 




( 27 ) 
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Figure 3. Gap Heat Transfer Coefficient 



The heat conducted out of the fuel is governed by Fourier's 
equation and is equal to the heat transferred across the gap 



’ = - k F3T < 28 > 

Equating equations (27) and (28) gives the fuel-clad inter- 
face condition 



k r- ( > z ) 3T p 

Tp ( a , z , t ) +p fij TF“( a ,z,t) = T c ( a > z > t ) (29) 



and from continuity 



3T p 3T 

kp ( a , z ) yp-( a , z , t ) = k c (a,z) -j^(a,z,t) 



(29a) 



As is common practice in convection heat transfer 
analysis, a coolant surface heat transfer coefficient, hsurf, 
is used to account for the thermal resistance at the clad- 
coolant interface. Similar to the fuel-clad interface analy- 
sis, the clad-coolant interface conditions may be determined 

k (b,z) 3T 

T ( b , z , t ) + yp —£(b,z,t) = T ( b , z , t) (30) 

c "surf 3r co 



and 



9T r 3T co 

k c (b,z) 7r^( b ’ z ’ t ) = k co (b ’ z) 7f^ (b ’ z>t) 



(30a) 



5 . Boundary Conditions 

The boundary and initial conditions for the heat 
transport problem are: 
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1) 



0 



37“ (°> z >t) 



2) 


T 

CO 


(r,- |,t) 


3) 


3T 

CO 
9 Z 


(r, j,t) 


4) 


8T co 

9r 


(c,z,t) 


5) 


8T f 

9z 


(r,+ !j-,t) 


6) 




T f (r,0) 


7) 




T c (1,0) 



8 ) T co ( ^ 0) 



t plenum 

0 



3T r H 

ST- 



= 0 



Tf(r) 

T°(D 



T° (r) 
c o v — ' 



Boundary condition 1) results from the assumed azimuthal sym- 
metry. The coolant has been assumed to enter the flow 
channel at a constant temperature, "'’plenum* This results in 
condition 2). Boundary conditions 3) and 5) result from an 
assumption that no heat is transferred axially from the fuel 
rod. Boundary condition 4) is the result of the use of the 
equivalent cell. Conditions 6), 7), and 8) are the assumed 
initial conditions. 

Solution of the nonlinear coupled neutronic and 
energy transport problem involves the solution of the partial 
differential equations (11), (12), (21), (24), (25), and 
(26), with the appropriate boundary, interface, and initial 
condi ti ons . 

Several works. Refs. [4], [11], [12], have demon- 
strated the feasibility and success of the finite element 
method in solving nuclear reactor dynamics problems. The FEM 
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is used to reduce the partial differential equations developed 
in this analysis to a system of ordinary differential equa- 
tions. Integration of these ordinary differential equations 
(ODE) yields the solution. 



C 
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IV. FINITE ELEMENT FORMULATION 



In this section, the basic theory underlying the finite 
element method is formulated. Selection of the finite ele- 
ments and the shape functions for the elements are given. 

Some simple transformations which facilitate the integration 
necessary in the FEM are also presented. 

A. BASIC THEORY 

To obtain a numerical solution, the governing partial dif- 
ferential field equations are transformed into a system of 
ODE in finite dimensional vector space. This may be accom- 
plished in several manners such as the finite-difference 
method, the variational method, or the weighted residual 
method. In this work, the Galerkin method (a weighted re- 
sidual method) is utilized for the discretization of the 
spatial domain. The Galerkin procedure will be applied to 
each of the governing field equations. Any of these equa- 
tions may be considered to be in the following form 



where represents the unknown function, e.g., in equation 

(21), ip represents <j>p , i represents the operator for 
each individual equation, and f is a forcing function. In 
the finite element method the solution is approximated as 



It (I’t) ~ iMu.t) = f (r,t) 



( 31 ) 



N 




( 32 ) 
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where 



N - the number of degrees of freedom 
N.j - the element shape function 
^ - unknown coordinate function 
< > - matrix notation for a row vector 
{ } - matrix notation for a column vector 
<N.> - <N -j N 2 ... N i ... N n > 




The residualj R(r,t) , is a measure of the error in this 
finite element approximation. The residual may be considered 
as 

R = || - - f (33) 

The best solution for $ is one which "minimizes" this 
residual. Various "minimums" are obtained by the weighted 
residual method by setting 



J W i (r) R dV = 0 i = l ,2, . . .N (34) 

V W.(_r) - weighting functions 

With the Galerkin method, the weighting functions are the 
shape functions defining the approximation of equation (32) 
(i.e., W. = ). A noteworthy attribute of the Galerkin 

method is the opportunity of using an i ntegrati on-by-parts 
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of the terms involving the second order spatial derivatives. 

A lower order finite element may be used than would have 
been possible otherwise. Once the weighting functions have 
been chosen, the problem becomes 

f N ■ {-|t - f} dV = 0 i = 1 , 2 , . . . N (35) 

V 

The integration involved in equation (35) is carried out 
on the element level, taking advantage of the use of a "local" 
coordinate system. Once the integration is accomplished, 
the results are merged into a system using "global" coordi- 
nates. On the element level 



$ e = <N . > e { 4; , } e j = l , 2 , . . . , N (36) 

J J 

where the superscript e indicates the element level. Sub- 

p p 

stituting $ into equation (35) and noting . } is not 

J 

a function of the spatial domain yields 



/<N 1 > e {N j .} e dV e {ii» j .} e -J* [<N i > e l {N j } e -<N i > e f e ]dV e {^ j } e = 0 
V V 

(37) 

where i , j - 1 , 2 , . . . N e 

N e - number of degrees of freedom for 
an element 



The operator i will vary depending upon which governing 
equation is under consideration. 



41 



B. SHAPE FUNCTIONS 

The shape functions, N.. , are chosen to satisfy certain 
completeness and convergence criteria [13] and will depend 
upon the finite-element used for the spatial discretization. 

Many previous works. Refs. [4], [11], and [12], utilized 
linear triangular shaped elements to discretize the spatial 
domain. This element was the first element considered. 
However, because the width of the cladding is very thin, 
elements in the cladding region would have extremely large 
aspect ratios (ratio of base to height) unless an extremely 
large number of elements in the axial direction were used. 

A large number of elements becomes numerically untractable. 
Previous experience with triangular elements had shown that 
large aspect ratios yield inaccurate results. Zlamal, Ref. 

[14] , showed the error, e, when using triangular elements, 
is proportional to the square of the longest side, h, and 
inversely proportional to the sine of the smallest angle, y 

e a h^/siny 

A triangular element with a large aspect ratio necessarily 

i 

must have a small related angle which adversely affects the 
error in the FEM. Hopefully to alleviate the problem, an 
isoparametric, quadratic, rectangular element was selected. 
The aspect ratio would still be large, but experience. Ref. 

[15] , with the use of rectangular elements indicated that a 
large aspect ratio is not always a detrimental factor. 
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The shape functions for this element are well documented. 
Ref. [13]. Utilizing a "local" coordinate system (See Figure 
4), the shape functions may be written as 

Corner nodes i = 1,3, 5, 7 

N i - i(l+? o )O*n 0 )(V\-') < 38 > 

Midside nodes N. = 1 ) ( 1 +n Q ) i = 2,6 (38a) 

N i = l( 1+ 5 0 )(1 " n2) i=4 ’ 8 (38b) 



where 




nn. 



These normalized shape functions are shown in Figure 5. 

The local and global coordinates are related by the 
following 

r = <N i > e { r^ } e (39) 

and z = <N\j > e {z^ } e (39a) 

C. COORDINATE TRANSFORMATIONS 

When using a local coordinate system, some simple trans- 
formations facilitate the integrations requ.ired by equation 
(37). In cylindrical coordinates with azimuthal symmetry 

dV = 2irrdrdz (40) 



The derivative terms may be transformed by the following 



I 3N . ’ 

— 

3z 




( 41 ) 
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Figure 4. Element Transformation 
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CORNER NODE 




Figure 5. Normalized Shape Functions 
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where [J] ^ is the inverse of the 2x2 Jacobian matrix de- 
fined in Appendix A. As shown in Appendix A, this inverse 
can be easily shown to be (for this problem) 



[J] 



-1 




(42) 



where 



★ ★ 



J 11 


= 2/Ar , 


J 12 = 


0 


★ 




★ 




J 21 


= o 


J 22 = 


2/ AZ 


Ar 


- radial 


length of 


the element 


AZ 


- axial 


length of 


the element 



Elements of area transform as 



drdz = det[J] d£dn (43) 

For this particular problem, det[J] may be shown to be 

(Appendix A) 

det[J] = -J- (44) 

where A e - area of the element 

Elements of axial length become 

dz = j- (45) 

where L e - axial length of the element 
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Utilization of these transformations makes the integrations 
required in equation (37) amenable to integration by numeri- 
cal Gaussian quadrature. 
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V. APPLICATION OF FEM TO GOVERNING 
FIELD EQUATIONS 



In this chapter, equation (37) is applied to the govern- 
ing field equations previously derived. The element matrices 
for each of the operators are developed so that the discreti- 
zation of the spatial domain may be accomplished. The inte- 
grations required by the application of equation (37) are 
performed numerically using Gaussian quadrature. 

A. GAUSSIAN QUADRATURE 

Prior to the application of equation (37), it is appropri- 
ate to discuss briefly the procedure used for the numerical 
integration. The product Gaussian quadrature formula is [16] 



The values of the weights associated with each Gauss point 
are given in Ref. [16]. Equation (42) may be simplified 
somewhat by combining the summations and weights 




.1 -1 



m m 



gU,n)d n d£ « Z E w i w j 9(5 i ,n j ) .(46) 

i = 1 j = l 



where 



1^ - area integration 

g(£,n) - any function of £ and n 

W. • - weight associated with location 

1 > J -i n v* s 



i or j 



m 



number of Gauss sampling points 
in one-dimension 




(47) 



k=l 
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where k - ixj 

W k - U, x Wj 

For line integrations, the Gaussian quadrature formula in- 
volves only one summation 

1 m 

I, = f f(n)dn = £ W. f( ni ) (48) 

_1 i = l 

B. NEUTRONIC FIELD EQUATIONS 

The discretization of the spatial domain by the finite 
element method is accomplished by applying equation (37) to 
the governing field equations, using 

4> k = < N j > 6 j = 1 , 2 , . . . , 8 (49) 

k= F , c , co 



1 . Fuel Region 

The governing equation for the fuel region, equation 
(20) after applying equation (37) becomes 

2l, f f r"5t- rdrd z-2’r//N 1 { f ^(rOp ,/) 

r z r z 

a ^ cp p 

+ fz (D F 3T 1 ) +Z aF^F (1 " e)[k - +p ' b1n(T F /T F )] 

t 

+ XB^aF f e' X i (t-t ' ) (k°+p)c() F dt , +XC 0 e" Xt }rdrdz = 0 
o 

(50) 
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The second order terms in equation (50) may be re- 
duced to a first order term by application of Green's Theorem 
or equivalently i ntegrati on-by-parts (See Appendix B). 
Dividing through by 2tt and reducing the second order terms 
yields 



/ rN i D F dz + / rN i D F rf ir *J J 



9<f>, 



F dr 

z r 



N i 8<f> F 

7~ Jt~ 



rdrdz 



r z 



r r ( Op Op 

+ f f | D F^3r~ T r~ + Jz~ TT^ -N i E aFM 1 C k ° + P- bl n ( T F /F Fo^ 



t _ X 

- N i A B E a f f e " X ( > ( k ° + p ) <|> p dt ' -N . AC° e" X 1 J 

n ' 



0 e _;)lt \ rdrdz = 0 ( 51 ) 



From continuity and boundary conditions the line in- 
tegrals are zero. Now using the approximate functions of 
equation (49) and noting that Op} 6 is not a function of 
space, equation (51) may be written as 



7 rdrdz{i|;p} e +^Tlp[{N nr } e <N jjr > e +{N i}Z } e <N j}Z > e ]rdrdz{^p} e 

^/^£ a p(l-B) C k OT +p_bl n ( T F/ T p) -I { ^i >e<N j >e r drd z { ^F }e -^^aF ' } ( k ^+ p ) d t‘ 



//{N. }<N .>rdrdz{^ c } e - AC°e~^ t ff {N-Odrdz = 0 
rz 1 J F rz 1 



(52) 
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The integrations may not be easily carried out with 
a shift to local coordinates and with use of the previously 



derived transformati ons . Rearranging equation ( 52 ) and 
assuming the properties are constant for each time step gives 

°F £ ^i i*}<^H*^j ,£ >+ ^i ,n^22.*}< J 22*Nj jr) >]rdet[J]d^dri{^p} e 

1 1 e 
-2 aF 0-6)(k°+p) / / {N i ><N j > rdet[J]dCdn(^p} e 

1 1 T F e 

-Z ap (l-B)b / ,/ In yt {N i )<N j .>rdet[J]dCd n {t|; p } e 



-A6Z aF f(t) 



1 1 
/ / 

-1 -1 



{N^ }<Nj>rdet[J]dCdn{^ F ) e -AC 



■At 



1 

/ 

•1 



/ (N.j }rdet[J]d£dn 



+ -7 / / (N.}<N.>rdet[J]dCd n {^} e = 0 
-1 -1 1 J *■ 



( 53 ) 



Since the element chosen has eight degrees of freedom 
(nodal points), the discretized matrices which result from 
the integration of equation ( 53 ) will be 8x8 matrices and the 
forcing function will be an 8x1 vector at the element level. 
Defining the matrices as 



1 1 



.{ _{ C' N i,c J ll )<J 11 N j,c >+{N i,n J 22 )<J 22 N j,n >]rdet[J] d?dn 



A 



e *ti 



*- H12 ij^8x8 4 2 ^( N i ,c) k (Nj,c) k J ll 

k = l 



*2 



+(N i J h ) k (N j,n ) k °22 3 r k W , 



( 54 ) 
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where 



H k - W jWj 



and 



r k = S r i (N i ) k 
i = 1 



1 1 

f f {N i }<N j >rdet[J]d5d n = [H3.^] 8x8 




( N i>k( N j)k^k W k 



( 55 ) 



11 e 

/ / {N.}rdet[J]dCd n = , = i E 

-1-1 1 i axi 4 k=1 



< N i>k r k W k 



(56) 



1 1 

/ / ln(T F /T D ){N i }<N j .>rdet[J]d?d n = 

2 

0 m 

= T 2ln(T F /T r ) k (N i ) k (N j ) k r k H k . (57) 

k=l 

To carry out the summation of equation (57), the 
temperature Tp must be known; however, the temperature is 
exactly what is being sought. To alleviate this problem, a 
linearization is used. 

In the solution technique, the temperature is pre- 
dicted for the next time step. It is this temperature which 
is used for the determination of matrix H4. 

Equation (53) simplifies to 

{D p [H12. j ]-[S aF (l-6)(k°+p)+X6i: aF f(t)][H3 ij ]+S aF (l-6)b[H4 ij .]}4 F } e 

-XC°c -Xt {F i } e + -^rCH3 - j ] (^p) = 0 



5 2 



(58) 



The function f(t) is evaluated by summing the 
values at each time step using the trapezoid rule for numeri- 
cal integration. 

f(t) = e" Xt /V^'C k° + p]dt’ (59) 

o 



Def i ni ng 

ij[g(t)] - j h i {g(t 1 ->+ gftj.,)) 
g(t) = e Xt < k°+p) 

h ^ - time step taken 
9 o = 0 



The function may be expressed as 

S 

f(t) = Cg(t)] 

i-= 1 



(60) 



(61) 



S - number of time steps 
2 . Clad Region 

Following the same procedure as with the fuel equa- 
tion, the discretized form of equation (11) becomes 



{D c [H12 ij D ♦ - >3 fj ]{j, c } e = 0 (62) 



3 . Coolant Region 

The governing equation for the coolant region, equa- 
tion (12), may be discretized into the form 

(“coCH^j] ♦ £ aco [H3 (j ]>{* co > e + ^CH3„K* co } e - 0 (63) 
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Once the governing equations have been discretized 



at the element level, they are combined into a system of 
equations at the global level. On the global level the 
governi ng equati on f or neutron transport takes the form of 



where [H] , [P], and [p] represent the system matrices and 
n is the number of nodal points used in the discretization. 
There are, then, n simultaneous ordinary differential equa- 
tions used to describe the neutron transport problem. 

C. HEAT TRANSPORT FIELD EQUATIONS 

The spatial domain for the heat transport problem is 
discretized in the same manner as the domain for the neutronics 
problem. The same element matrices previously defined are 
valid. Let 



1 . Fuel Region 

The governing field equation for the fuel, equation 
(24), is discretized by applying equation (37). Using an 
integration by part to lower the order of the second order 
terms allows equation (24) to be written as 




T k -- <V e{T kj }e J" 1 - 2 -- 



8 



(65) 



k = F , c , co 



r 




3N. 9T 
3z“ Iz 



( 66 ) 



54 



From continuity considerations the line integrals 
are zero except along the boundaries of a region. It is 
assumed that no heat is transferred from the cell in the 



axial direction (boundary condition 5); therefore, the 
first line integral of equation (66) is zero. In the neu- 
tron diffusion problem there was continuity of flux at the 
interfaces so that the line integrals were zero; however, 
the heat transfer at the interfaces is affected by the gap 
and film conductances. The fuel-clad interface condition, 
equation (29), may be rewritten as 



Substituting equation (67) into (66), dividing by minus one 
and utilizing equation (65) yields 




c 



(67) 





r z 



/ / eE fF (N i } e <N j > e rdrdz{^ F } e 



+ / / PpC pF (N i } e <N j > e rdrdz {T p } e = 0 



r z 



(63) 



e* - element across the interface 
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Transforming to local coordinates and integrating 
by Gaussian quadrature yields the same e 1 emen t. ma tri ces as 
given for the neutron flux, equations (54) and (55), except 
for the line integrals between regions. It should be noted 
that the line integrals exist only on the interfaces, along 
which there is a discontinuity of temperature. For the 
fuel equation the interface corresponds to the local co- 
ordinate £ = 1 . Define 



1 

/ (N.}<N.>dn 

-1 1 J 



vn 



^ k1 ij^8x8 



? £ < N iV N A w k 

k=l 



(69) 



where the N's are evaluated at £=1. Many of the terms of 

t 

K1 will be zero since only the nodes on the C=1 boundary 
will have shape functions which are non-zero. It is through 
K1 that the temperatures for each region are coupled to- 
gether. Equation (68) may now be written as 



r a H gap {CK1]{T F }e ' [K1]{T c }e * } + k F [H12:|{T F }e 

-eE f p [H3] { ip F } 6 + P F C pF [H3]{x F } = 0 (70) 

In obtaining this equation, it was assumed that material pro 
perties for each nodal point were constant at each time step 
Perhaps a better assumption would have been to assume an 
average value for the properties of each element. The dif- 
ference should not be significant, and the assumed constant 
nodal properties were numerically more tractable. 
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2 . Clad Region 



Applying equation (37) to the governing field equa- 
tion for the clad, equation (25) gives the discretized form 
of the equation. Assuming no heat transfer in the axial 
direction on the boundaries (Boundary condition 5), equation 
(25) becomes 



3T 



- ! [N i rk c aria dz V 1 rdrdz 



3N. 3T 3N_. 3T 



r z 



+ / / P r 
r z c 




3T 



i 3t 



rdrdz = 0 



(71) 



In the cladding, there are two interfaces along which 
the line integral of equation (71) is not zero, along the 
fuel-clad interface and along the clad-coolant interface. 

For the fuel-clad interface, equation (61) applies. For the 
clad-coolant interface, the interface condition, equation 
(30), may be rewritten as 

k c 5T ' hi " rf (V T c.l ' -‘co r 5 < 7Z > 



Along the fuel-clad interface, the local coordinate 
corresponds to £ = -l. Define the new element matrix 



/ / {N.} e <N,> e dn 

- 1-1 1 J 



, e m 

^8x8 "7 E < N i>k <Vk W k 



k=l 



where the N's are evaluated at 5=-l. As with K 1 , K2 will 
have many zero values because the shape functions are 
eval uated at £ = - 1 . 
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Along the clad-coolant interface, the local coordi- 
nate corresponds to £=1 and the K1 matrix is appropriate. 

Substituting equations (67) and (72), equation (71) 

becomes 



/KurfVWIb dz -/K,ap N i (T F- T c>la ^ 



3N. 3T 3N. 3T 

+ 1 ! k c [ 3F IT + ST 3F- ] rdrdz 



r z 



+ ' / d c C pc "l TT rdrdz s 0 



(73) 



The governing equation for the clad region may now be written 
as 



+ r b h surf { [ K1 3 (z c )e -[KlKT co > e *> 

+ k c [H12]{T F } e + p c C pc [H3]« c ) = 0 (74) 

The line integrals of equation (73) affect only nodes which 
are on one of the boundaries; therefore, the nodal inputs 
into K1 and K2 are zero unless the node is on one of the 
boundari es . 

3 . Coolant Region 

The field equation governing the coolant may be 
discretized in the same manner as above. After applying 
the Galerkin method and performing an integration by parts 
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on the second order terms, equation (25) 
becomes 



^ rN i k co &l ' dz * // k co l 3r 



3N. 3T 3N . 3T 
r i co_,_ i co i 

L*- 3r + 3z 3z -* rdrdz 



+ 



J ^ V co p co C pco N i 3z 



co 



rdrdz 



+ 



/ / P 
r z 



co 




3T 
N i 3t 



co 



rdrdz = 0 



(75) 



The line integral, when evaluated at c, is zero (boundary 
condition 4). When evaluated at b, or correspondingly at 
£ = -l , equation (72) is valid and K2 matrix is appropriate. 
All the terms of equation (75) have been defined except the 
flow term. Define 

_/ 1 _/ 1 {N i } e< N ijn >er det[J]d£dn = [H5, j ] 0x8 

A e m2 

= ? £ (N i } k (N i,n } k r k w k (76) 

k = 1 



Transforming to local coordinates and integrating reduces 
equation (75) to 



r b h s u rf<[ K2 ^co }e -[ K2 :i {T c >e > + k co [H,2]lT co }i 



♦ VcoSco^co 1 * + p co C pco [H3]{i co >e ' 0 < 77 > 
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Now that the governing equations have been defined 
for each region on the element level, equations (70), (74), 



and (77), they may be assembled into a system equation on 
the global level. The equation will be in the general form 
of 

^^nxn^nxl + ^^nxn^^nxl + ^ ^nxn^nxl " 0 

(73) 



D. DISCRETIZATION OF THE SPATIAL DOMAIN 

Prior to the numerical solution of the governing equa- 
tions, equations (64) and (78), the spatial domain must be 
divided into a number of elements. For this work the domain 
was subdivided as shown in figure 5. 

Since there is a discontinuity of temperatures at the 
interfaces, as described by equations (29) and (30), a 
noval application of the FEM method was necessary. The com- 
mon practice for handling these ‘'flux" type boundary condi- 
tions is to define a constant reference temperature, T^ , 
as when working with a convection heat transfer problem [17], 
or to define a known function, as when working with a 
fracture mechanics problem [18]. In either case the refer- 
ence condition was known. The novel application here lies 
in the use of a different field equation to describe the 
reference temperature, e.g., the clad equation (74) is the 
reference condition for heat transfer from the fuel across 
the gap interface. 
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The discontinuity of temperature at the interface neces- 
sitated another novel application of the FEM. Since there 
is a temperature drop along each interface, a single node 
there is not adequate. In the discretization of the domain, 
two nodes were used for each interface point (for example, 
points 62 and 63 in figure 6). This allows the temperature 
drop due to the gap and film conductances to be taken into 
consideration. Since two nodes are used, the governing equa- 
tions for each region are not directly coupled together. 

The coupling of the regions is accomplished by the "flux" 
boundary or interface conditions since it is assumed that 
any heat flux leaving a region enters the adjacent region. 
Consider a typical set of elements on an interface 
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FUEL 




Figure 6. Finite Element Discretization 
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The coupling terms K1 and K2 may be combined into a system 

K matrix which shows the coupling. The K matrix for the 

simple set shown is 

1 2 3 4 5 6 7 8 • • • 13 14 15 16_ 

a = ’6/15 
b = 4 / 15 

c = V 15 

d = V 15 



As can be seen, the nodes on the interface are coupled to 
the adjacent element interface nodes. For example, node 2 
in element I is coupled to nodes 3, 8, and 14 in element II. 

E. OPTIMUM COMPACTING SCHEME 

The system matrices ( H, P , etc.) are nxn matrices, 
where n is the number of nodal points used in the discre- 
tization of the domain. In terms of computer storage, these 
matrices may become excessively large if they are stored as 
nxn. There are several techniques available to reduce this 
storage requirement. The most common method is the banded 
storage scheme, whereby only the banded portion of the 
matrices are stored. With judicious numbering of the nodes. 



1 

2 

3 

4 

5 

6 
7 



a -a 
- a a 



C -c 
c c 



c -c 
-c c 



b -b 
-b b 



C — € 
-C C 



d -d 
d d 



13 

U 

15 

16 



C -c 
C € 



d -d 
d d 



b -b 
b b 
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considerable savings may be realized. However, it is not 
the optimum storage scheme [8], 

J_ L. 

Since the shape functions, N i , for the k— nodal equa- 
tion are nonzero over only the element containing k, the 
system matrices are not only banded but sparse as well. 

The sparseness is due to the non-consecuti ve numbering of 

t h 

the nodes surrounding the k— node. The optimum compact 
storage (OCS) scheme compacts the matrices by storing only 
the non-zero elements of the matrices. The implementation 
of the OCS scheme requires two additional integer arrays, 
say JA and NAME. The NAME array identifies the nodal points 
which contribute to each nodal equation. The JA array acts 
as a pointer to indicate where the nodal equation starts in 
NAME. Consider the following simple 2x2 system with nodes 
as indicated 



7 8 

III 

r 4 5 


9 

IV 

6 


1 

1 2 


n 

3 



The NAME array starts with node 1 and identifies the nodes 
which contribute to node 1 (i.e., 5, 4, and 2). The NAME 
array would, then, give the nodes contributing to node 2 
and so forth, so that 

1234 • 5673910 *11 1213 14 ♦ • JC 

<NAME> = <1542 *. 254163 13 6 5 2 : ... .* 9856> 
and i 2 3 

<JA> = <1 5 11 ...> 
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The algorithm to assemble the element matrices into a compact 
storage vector is strai ghtforward and represents a significant 
savings in computer storage [8]. The system matrices are 
stored as a vector rather than a two-dimensional array. For 
example, the value which would be stored in position (1,5) 
of the nxn array is stored in position 2 of the system vector. 
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VI. NUMERICAL SOLUTION 



This section contains a brief description of possible 
solution techniques in addition to the solution technique 
chosen. Computer subroutines necessary to implement the 
technique are also described. 

A. SELECTION OF METHOD 

The numerical solution of the system of implicit ordi- 
nary differential equations, equations (64) and (78), may 
be accomplished by any of a number of different techniques 
such as Houbolt's method, Crank-Ni colson' s method. Gear's 
method, or implicit Gear's method. It was not the objective 
of this analysis to determine which of the numerical solu- 
tion schemes is the most efficient. Each method has its 
advantages and disadvantages. The Cra nk-Ni col son method 
is a single-step, implicit equation solver and, therefore, 
does not require storage of previous time solutions. When 
analyzing neutronic problems, the system of equations which 
arises is commonly very stiff (i.e., a rapid change in flux 
over a short period of time). The Crank-Ni col son method 
has, in a past work [8], demonstrated difficulty in tracking 
these stiff systems. Gear's method was specifically 
developed for stiff systems and can handle the problem very 
well. However, Gear's method is a multi-step, predictor- 
corrector method requiring storage of previous time solutions. 
In addition to this disadvantage. Gear's method requires the 
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transformation of the developed implicit O.D.E.'s into an 
explicit system of O.D.E.'s. After this transformation is 
done, the system matrices are, no longer sparse or banded, 
thus eliminating the use of the optimum compacting scheme. 

In an effort to overcome these difficulties. Gear's method 
was modified. Ref. [9], to treat the implicit system of 
equations as well as to allow use of the optimum compacting 
scheme. A previous work, Ref. [8], has shown that the 
implicit Gear's method is particularly attractive in solv- 
ing the type problem developed in this analysis. Therefore, 
the implicit Gear's method is used for the solution of the 
system of O.D.E.'s arising i n thi s analys is . 

No attempt will be made here to give the mathematics in- 
volved in developing the implicit Gear's method. Reference 
[9] may be consulted if details are desired. A listingof 
the computer program developed will be given in the Computer 
Program section. In order to utilize the implicit Gear's 
method, several user supplied subroutines must be developed: 
1) DIFFUN, 2) JACMAT, and 3) NUITSL. 

B. USER SUPPLIED SUBROUTINES TO IMPLEMENT 

THE IMPLICIT GEAR'S METHOD 

1 . DIFFUN 

Subroutine DIFFUN evaluates equations (64) and (78) 
for a given time and for given values of ip , $, x and x . 
Since at each nodal point, i, there is a solution for the 
flux and for the temperature, the solution was set equal to 
DYI and DY II, respectively. In addition to having flux and 
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temperature at each nodal point, there are also three differ- 
ent regions in the domain which have different governing 
equations. An integer array, ITYPE, was developed to indi- 
cate for each nodal point whether it was: 0) a fuel node 

not in an interface element, 1) a fuel node in an interface 
element, 2) a cladding node or, 3) a coolant node. Using 
ITYPE, the computer program is directed to a different sec- 
tion depending upon the type of node being considered. After 
all the nodes have been considered, boundary conditions are 
established by changing DYI and DYII for the appropriate 
boundary nodes. Since in this analysis there is continuity 
of flux at the interfaces, special considerations must be 
given to these nodes. At the fuel-clad interface, the value 
of DYI for the clad node was set to the value of the flux 
at that node minus the value of the flux at the adjacent 
node (i.e., DYI.. = ip^ - i|>. ^). Similarly, at the clad-cool- 
ant interface, the value of DYI for the coolant node was 
set to the value of the flux at that node minus the value 
of the flux at the adjacent node. During the solution of 
the problem, DYI is driven toward zero, which in the limit 
forces ip . to equal ’J'-j.] • This is the desired continuity 
res ul t . 

2. JACMAT 

Subroutine JACMAT, evaluates the Jacobian matrix 
(for Gear's method) at the given time and for the current 
values of the dependent variables. The Jacobian for an equa- 
tion of the type, 

F(y, y, t) = 0 (79) 
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may be represented as [19], 



J 



r 3 F a 0 3F-, 



(80) 



where a Q and g Q are coefficients from Gear's method and 
h is the time step. Using the notation of DIFFUN, let DYI 
and DYII represent equations (64) and (78), respectively. 
The Jacobian matrix may, then, be written as (J is called 
PW in JACMAT.) 



PW 



r 3DYI a o 3 DY I 3DYII a o 3DYII-, 
L9 * '67 si , ’ 9t ’ 3 . J 



(81) 



It is the form of equation (81) which is programmed in 
JACMAT. As in DIFFUN, the integer array ITYPE is used to 
indicate the appropriate section of the program to be 
utilized. The problem boundary conditions must also be 
accounted for in JACMAT. In DIFFUN, the value of DYI or 
DYII was set to zero for constant boundary conditions (i.e., 
zero). This cannot be done in JACMAT since a division by 
zero would occur. For a constant boundary condition at the 
ith node, the value of PW is set to one for the diagonal 
term and zero for all other terms of the ith equation. 

3. NUITSL 

Subroutine NUITSL solves the system of equations 
for the quasi-Newton iterates. In this analysis the sys tern 
is solved using a successive over-relaxation (SOR) method. 

In this work, the optimum amount of over-rel axati on was not 
determined. Since no effort was made to find the optimum. 
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it was felt a small over-relaxation would be best. The over- 
relaxation factor of 0.02 was used. For small values of 
this factor, the SOR method approaches the Gauss-Siedel 
iteration technique. 
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VII. PROCEDURE 



In this section, the method utilized to obtain a solu- 
tion is described. The input data necessary to run the 
developed computer program will be documented. 

Prior to initiating a transient overpower excursion, the 
steady-state conditions for the fuel cell must be known. 
Since thesystemof equations which were developed are not 
specifically designed to obtain a steady-state solution, the 
initial steady-state conditions must be part of the input 
data. The initial temperature distribution was obtained 
from the steady -state conditions given in Ref. [7 ] . 

The axial temperature distribution for the fuel center- 
line, fuel surface, clad, and coolant have been determined 
for several different fuel life cycles [7], For this analy- 
sis, the beginning of life cycle for channel 10 was used. 
Although this distribution is somewhat artificial, it should 
be adequate for this analysis. It is the trends of the 
results which are considered important. The distribution 
within the fuel radially is taken to vary as the square of 
the radial distance; then 

2 2 

fp( r J z jO) = Tp(0,z,0) (1- + T p ( a , z , 0 ) 

d d. 

Within the cladding and the coolant, the initial radial tem- 
perature distribution is assumed to be constant. 
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The initial flux distribution is assumed to be radially 
constant, a flat flux assumption. In the axial di recti on 
the flux is assumed to vary as the shape of the sine func- 
tion. The maximum flux, the flux at the axial center, is 

an input parameter. For this analysis, the maximum initial 

1 4 2 

flux was taken to be 10 neutron/cm sec. 

To obtain a steady state flux distribution, the value 

of fission cross section for the fuel is varied. A trial- 

and-error method is used until a critical fission cross sec- 
c r 

tion, , which gives a steady flux is obtained. 

Once the steady-state conditions have been determined, 
the excess reactivity may be inserted. This starts the 
transient overpower excursion. 

A. INPUT DATA 

The first data card contains: the order of Gauss quadra- 
ture, the number of radial elements in the fuel, the number 
of axial elements, number of nodal points in the radial 
direction, and the height of the fuel rod. The next cards, 
one for each radial nodal point, contain the nodal radial 
distances. The next cards contain the fuel centerline, fuel 
surface, clad and coolant temperatures. There is one card 
for each axial node. The next card contains the maximum 
flux. The next four cards contain the physical parameters 
listed in Table I. The last input cards contain the time, 
end time, estimated initial time step, minimum time step, 
and maximum time step. A sample data deck is shown in figure 
7 . 
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TABLE I. Physical Parameters 



Fuel Diffusion Coefficient (DCF) 

Doppler constant (B) 

Energy release per fission (E) 

Fraction of delayed neutrons (BETA) 

Fraction of delayed neutrons for 
the ith group (BETAI) 

Decay constant for ith delayed neutron 
group (DCLAMI) 

Initial flux for delayed neutrons 

Average number of neutrons released 
per fission (ANU) 

Neutron velocity (VEL) 

Fuel absorption cross section (SIGAF) 

Critical fuel fission cross section 
(SIGFF) 

Blanket absorption cross section (SIGAB) 
Blanket fission cross section (SIGFB) 

Step reactivity input (RHOA) 

Ramp reactivity input (RHOB) 

Fuel density (DENF) 

Clad diffusion coefficient (DCC) 

Clad specific heat (CPC) 

Clad density (DENC) 

Clad thermal conductivity (TKC) 

Clad absorption cross section (SIGAC) 
Coolant diffusion coefficient (DCCO) 
Coolant absorption cross section (SIGACO) 
Coolant flow velocity (VCO) 

Surface heat transfer coefficient(HSURF) 



0.93 [cm) 

0.006 

7.652xl0~ 13 [cal/fissions] 
0.0064 

0.0064 

0.0784 

10 2 
1x10 [neutron/cm sec] 

2.44 

4.8 xlO^ [cm/sec] 

0.088 [cm" 1 ] 

0.0586875 [cm -1 ] 

0.0 [cm -1 ] 

0.0 [cm" 1 ] 

Variable 
Variable 
10.9 [gm/cm 3 ] 

1.1 [cm] 

0.12 [cal/gm °C] 

8.0 [gm/cm 3 ] 

0.0526 [cal/cm sec °C] 
0.0015 [cm' 1 ] 

1.55 [cm] 

0.00004 [cm' 1 ] 

396.0 [cm/sec] 

0.7 [cal/cm 3 sec °C] 
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VIII. RESULTS 



When using the finite element method, one of the first 
considerations must be given to the convergence of the 
method. To determine convergence, the results for a given 
point are compared for different finite element discretiza- 
tions. As shown in figure 8, the results are comparable 
but no definite claim of convergence can be made. However, 
for this work, it was felt that these results were adequate. 

It was not the object of this analysis to arrive at the 
"final" result; it was the trends and methods that were of 
interest. Since the 66-element mesh appears to give a fair 
approximation of the results, the 66-element mesh was used 
as the discretized domain. 

The next item of consideration was the determination of 
a neutronic steady-state condition. This proved to be a 
very time consuming task. The fission cross section for the 
fuel was varied by a tri al -and-error method in an attempt 
to find the critical cross section which would give a steady 
state. As may be seen in figure 9, a change in cross sec- 
tion of less than one percent significantly affected the 
state of the problem. It was felt that the critical value 
was between the values of 0.05875 and 0.058625. Time did 
not permit investigation for critical value; therefore, 
it was assumed that the value for the critical fission cross 
section was half way between the values (i.e., z£ r = 0.0586875) 



75 




76 



Figure 8. Convergence of Finite Element Method 
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Fi gure 



Even if this value is in error, which it most probably is, 
the net effect would only be a small decrease or increase 
in the proposed reactivity insertion. The reactivity inser- 
tion would, then, be an approximation of the actual reactiv- 
ity of the problem. 

The first test problem considered was a step increase 
in reactivity of approximately ten dollars. For the uranium 
dioxide fuel, one dollar of reactivity was taken to be 
0.0064. Figure 10 shows the time history of the flux at 
the center of the fuel rod. Figure 11 gives the correspond- 
ing temperature profile. The temperatures were taken at the 
hottest point of each region at the axial center (i.e., the 
fuel centerline, clad inside surface, and coolant inside 
surface). As seen on the fuel temperature time history, the 
fuel rapidly reaches the fuel melting point. The model 
developed does not take into consideration melting of the 
fuel. This melting would tend to decrease the effect of the 
transient. The problem was allowed to continue despite this 
inconsistency in the mathematical model. A short time after 
fuel melting, the inside surface of the clad reaches its 
melting point. The temperature in the coolant experienced 
what is felt to be a numerical phenomenon. The coolant tem- 
perature decreased prior to the small rise at the end of the 
transient. Intuitively, this decrease does not seem to be 
realistic. A similar occurrence was observed while conduct- 
ing sample tests on the developed computer program. Refer- 
ence [20] reported the same phenomenon. It is felt that this 
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Figure 10. Flux Profile 
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Temperature Profile 



phenomenon is a quirk of the finite element method. The 

radial and axial distributions of the neutron flux are pre- 

_ 3 

sented in figures 12 and 13 for time equal to 7.39 x 10 
seconds. The distributions are, basically, as anticipated. 
The neutron flux peaks slightly before the axial center. 

It was expected to peak at the axial center. The radial 
and axial temperature distributions for the same time are 
presented in figures 14 and 15. As with the flux, the tem- 
perature profiles were, basically, as expected. The fuel 
temperature peaks slightly below the expected location, 
most likely in response to the peak in axial fluxes. The 
coolant unexpectedly drops near the outlet of the fuel rod. 
The finite element method characteri s ti cal ly has some prob- 
lems on the boundaries of the domain; this may account for 
the drop in coolant temperature. 

The temperature of the fuel and cladding do not appear 
to be as closely coupled as anticipated. As seen in figure 
14, a significant increase in fuel temperature has resulted 
in a relatively small clad temperature increase. As noted 
in figure 11, there appears to be a time lag in temperature 
response for each region which may account for part of the 
apparent temperature disagreement. It is felt that the 
temperatures should be more closely related, which indicates 
a higher gap heat transfer coefficient should be utilized. 
Since values of the gap heat transfer coefficient were 
assumed, it is not unreasonable to believe the values used 
are too low. 
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FIGURE 12 



RADIAL FLUX 
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Figure 12. Radial Flux 
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Figure 13. Axial Flux Profile 
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Figure 14. Radial Temperature Profile 
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Axial Temperature Profile 



Time did not permit investigation of other reactivity 
insertions. Other reactivity inputs may be investigated by 
students in the future. 

This work does not represent a solution to the very com- 
plicated nuclear reactor problem. It does represent an 
application of a numerical technique which is relatively new 
to nuclear applications. Methods for implementing the finite 
element method have been discussed, and a computer code has 
been developed for the simplistic model considered. 
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IX. RECOMMENDATIONS 



For the model developed, perhaps the most important item 
to pursue is the critical fission cross section. A better 
determination of this value is necessary so that the reactiv- 
ity insertion is more accurately known. Different test cases 
for the prompt critical and prompt subcritical reactor could 
then be conducted. 

In further developing the model, more consideration 
should be given the gap heat transfer coefficient. As noted 
in the results, the value used appears to be too small. 

Sample problems for different gap heat transfer coefficients 
would give a better indication of the values to use. 

Melting of the fuel during the transient would probably 
be the next major improvement on the model. With relatively 
few changes, the model could be adapted to allow melting 
element by element. This, too, would be an approximation 
but, still, an improvement to the model. Perhaps at the 
same time, a simplified model to take into consideration the 
fuel res tructuri ng could be implemented. 

Another improvement would be to consider reactivity 
feedbacks in addition to the Doppler feedback. Sodium void- 
ing and fuel rod expansion are two of the more important 
feedback effects to consider. 

On the numerical side, probably the most important thing 
to do would be to run the computer program on the 
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"H-compiler", which optimizes the program. However, on 
several runs using the H-compiler, erroneous results were 
obtained. With sufficient time, this could be corrected 
to allow use of the H-compiler. The use of the H-compiler 
results in a savings in computer time. The present program 
runs on a "6-compiler" and takes excessive amounts of com- 
puter time (two to four hours per run). 

In addition to this, the optimum over-relaxation factor 
in the implicit Gear's method could be determined by trial- 
and-error . 

Implementation of these recommendations should enhance 
the analysis and lead to a more efficient computer code. 



88 



APPENDIX A 



DEVELOPMENT OF TRANSFORMATIONS 



The Jacobian matrix [J] may be written (8) for two 
dimensions as 

N N 

E E N i,c z i 



[J] 



i = 1 
N 

I 

i = 1 



N. r. N. z. 

r-J i » n i " i.oi 



i = l 
N 

i =1 



(Al) 



For a simple 2x2 matrix [A] the inverse is 



[A] 



[A] 



a b 

c d 



-1 1 



de"t[Aj 



] 



-b 

a 



Applying this fact to equation (Al ) gives 

N N 



1 1 



[J] = de'tCJj 



^ N i,c z i 



i = 1 
N 





Eh, 

i = 1 1 


★ 




1 1 


J 12 


* 


★ 

"J 


21 


J 22 



i = 1 
N 



2 N. r- 
i = 1 1 ’ ^ 1 



( A2 ) 
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From matrix algebra 



det[J] 



= 2X 

i = 1 



C r i 



N 

E 

i =1 



z . 
i n i 



N N 

E n 1 i 5 z 1 

i = 1 i = 1 



(A3) 



The derivatives of the shape 
equati ons ( 38) . 

N i ^ = j(l+n) (2?+n) 

N 2,C = -^ 1+ n) 

n 3,c = i( 1+ n)(2s-n) 

N 4,c = - T* 1 '" 2 ) 

N 5,c = 7^ 1 " n ^ ( 2 ^ +ri ) 

N 6,C = 

N 7,C = iO-n)(2C-n) 



functions may be found from 

N ljn = j(l+C)(2n+C) 

N 2 ,n 3 I<>-5 2 ) 

N 3,n - }0-C)(2n-C) 

N 4 _ n - - n(l-C) (A4) 

N 5,n ' ?<>-5)<2nH) 

N 6,n * - 

n 7;I1 - }(l+C)(2n-5) 

N 8,n - "'( 1+ 5) 



If one now considers an arbitrary element 

z I 

< k 




with the midside nodes exactly at the midpoint (not a neces- 
sary criteria for the FEM) and substitutes into equation (A3), 
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the result is 





det[J] = 


( z 3 - z 1 > ( r 3 ‘ 
4 


• r i) . 


A e 

4 


(A 5 ) 


Substituting 


into equation (A 2 ) and 


using 


( A 5 ) will 


yield 




★ 

J n * 


A e [ 2 <Z3 ‘ : 


■ 


2/r 3 -r 1 , 


(A6) 




★ 

J 12 * 


0 , 






(A 7 ) 




★ 

J 21 ■ 


0 , 






(A8) 


and 


★ 

J 22 


“e 


-,)] ■ 


2 / z 3 ~ 2 1 • 


(A 9 ) 

/ 



The inverse of the Jacobian matrix now becomes 



[J] 



-1 



2/r 3 -r 1 0 

0 2 / 2 3 ~ 2 -]- 



(A10) 



For integration along a line, the transform at ion used 



1 s 



dz = det[J']dn 



(All) 



In this case 



det [ J ' ] = N • z ■ 
i = 1 



( A 1 2 ) 



Again considering the arbitrary element and substituting 
(A 4 ) into ( A 1 2 ) will result in 



det[J ' ] = 



z 3 ~ 2 1 L e 



( A 1 3 ) 
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APPENDIX B 



REDUCTION OF SECOND ORDER TERM 



The second order term of the governing field equations 
may be reduced to first order by integration by parts. 
Consider, for example, 

// N 4fF< rD !f> + f? (D i>] rdrdz (B1 > 

r z 



which may be expanded as 



If 



r z 



+ V If i + ",® f + $ * V 1? If] ««* 



(B2) 



Integrating just the second order terms by parts will yield 



If f 



D i^ rdrdz =j [^.rD-g^dz 
r z 9r z 



9N 

• / f lfC DN i + rN i I? * r[) IT - ] drdz < B3 > 

r z 



and 



r r 2 9N 

JJ N i D ff '' drdz -f [ N i D Hlz rdr D+N i $ rdrdz < B4 > 

r z 32 r r z 



Substituting the results of (B3) and (B4) into equation (B2) 
will give 



/[NirD fjf| r dz ♦/ [N.D |f| z rdr -fj D[^l $ * ^ |§] rdrdz (B5) 



r z 
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APPENDIX C 



LIST OF RELATIONS FOR MATERIAL THERMAL PROPERTIES 

A. FUEL (U0 2 ) 

1. Specific Heat, Ref. [19] 

C pF = [18.45+2.431xl0' 3 T-2.272xl0~ 5 T 2 ]/270.07 

[cal/gm °C] 

T - °C 

2. Thermal Conductivity, Ref. [5] 

Tk F = [l-2.5(l-p TD )]x[ 1 3 3 ;j + 4.79 x10" 13 T 3 ]x0.239 

[cal/cm sec °C] 

P TD “ P ercent theoretical density 
T - ° K 



B. CLAD (Stainless Steel) 

Properties are assumed to be temperature independent, 
and average values from Ref. [5] were used for the 
clad properties. 



C. COOLANT (Liquid Sodium) 

1. Specific Heat, Ref. [5] 

C = 0. 34574-0. 79226x10" 4 T+0.34086x10" 7 T 2 
P co 

[cal/gm °C] 

T - ° F 

2. Density, Ref. [5] 

p co = [59. 566-7. 9504x1 0" 3 T-0. 2872x1 0" 6 T 2 

+ 0.06035xl0"^T 3 ] x 0.01601 [gm/cm 3 ] 

T - ° F 
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3. Thermal Conductivity, Ref. [5] 

Tk = [54.306-1 .878 x10“ 2 T+2.0914x10" 6 T 2 ] 
co J 

x 4.134xl0~ 3 [cal/cm sec °C] 



T - °F 



D. SURFACE HEAT-TRANSFER COEFFICIENT 



'surf 



T k De V co p co C p n o 

co [7.0+0.025 ( — ) 0 • 8 j 



De 



co 



[BTU/hr ft^ °F] 



Tk co - [BTU/hr ft °F] 

p co “ [lbm/ft 3 ] 

V CQ - [ft/sec] 

C - [BTU/lbm °F] 
p co 

De - equivalent diameter [ft] 
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S1CRE YIT — THE INITIAL TEMPERATURE DISTRIEUTION 
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SUBROUTINE IN I T < I EP , N ELCON , R, Z , Y ) 

4* INITIAL SETS THE STARTING VALUES 
** FCF THE TEMPERATURE ANC THE NEUTRCN FLUX 
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TFE MAXIMUM FLUX (I.E. TFE FLUX 
AXIAL CENTER) RAOIAL FLUX WILL BE ASSUMED 
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SIERCUTJNE SHAPE 
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SIEFOUT INE MATK1 ( NELCCN ) 
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SLBRCUT INE JACOB ( R , N ELCCN , P , Z ) 
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SLERCUT INE MATH 12 ( y ,N ELCGN , R , Z ) 
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S l E RCUT I NE MATH4 ( N , N ELCQN , R , l , Y 1 T ,Y ) 

SLEFCUTINE MATH4 CALCULATES THE 8X8, TEMPERATURE 
CEFENOENT ELEMENT H4 MATRIX 
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